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STEADY INVISCID TRANSONIC FLOWS OVER PLANAR AIRFOILS - 
A SEARCH FOR A SIMPLIFIED PROCEDURE 

By R. Magnus and H. Yoshihara 
General Dynamic s/convair Aerospace Division 
San Diego, California 


SUMMARY 


A procedure to calculate steady inviscid flows over a general class of planar 
airfoils with embedded shocks at high subsonic Mach numbers has been developed 
earlier using an unsteady finite difference analogue. Here the embedded 
shocks were assumed to be sufficiently weak to permit the assumption of isen- 
tropic flow. The desired steady flow was obtained as the asymptotic flow 
for large times by prescribing suitable boundary and initial conditions. The 
unsteady Euler equations have been used in the proper conservation form to 
obtain the proper jump conditions across the embedded shock following the 
guidelines set forth by Lax using the concept of a weak solution to a hyper- 
bolic initial value problem. The Lax-Wendroff two step second order differ- 
ence scheme has been employed. The procedure as described above yielded 
satisfactory solutions for even the extreme cases of a highly blunt nosed 
aft cambered supercritical profiles, but at a cost of large computer times. 

In the present study simplifications are investigated to reduce the computer 
times without sacrifice of the accuracy of the solutions and always conform- 
ing to Lax's guidelines. As a first step the adequacy of the Guderley- 
Von Karman small disturbance theory with fictitious times derivatives and 
various simplified conditions at the airfoil have been examined to treat 
blunt nosed aft-cambered supercritical profiles. The results indicate that 
the simplifications are simply too severe resulting in an inacceptably poor 
prediction of the flow about the nose and. the shock wave. A more exact 
formulation of the problem is then proposed where the flow domain is con- 
sidered in two parts, the first sub-domain embedding the shock wave, while 
the second containing the remaining region where the flow is continuous. In 
the latter region the exact equations in non-conservation form for the two 
cartesian velocities are used, while the proper conservation form is employed 
in the sub-domain containing the shock. The exact flow tangency condition 
is prescribed on a semi-infinite slab where the upstream extremity is a semi- 
circle which overlays the airfoil leading edge. Downstream of the leading 
edge the condition is fulfilled in a quasi-planar fashion. The procedure is 



then illustrated by several examples with a supercritical profile. The 
report is then concluded by assessing the treatment of the embedded shocks 
both by the present unsteady and other steady finite difference analogues, 
suggesting here the need to incorporate viscous effects. 

1. INTRODUCTION 

In an early effort (Refs. 1-3) a finite difference procedure was developed 
to treat high subsonic flows with shocks over a general class of planar 
airfoils. Here the desired steady flow was obtained as the asymptotic flow 
for large times in an unsteady formulation for which a marching process 
could be reliably employed. The finite difference analogue here was based 
upon the unsteady Euler equations assuming the shocks to be sufficiently 
weak as to permit the assumption of isentropic flow. The airfoil condition 
was fulfilled exactly at the actual airfoil surface. Following the results 
of Lax (Ref. 4) the equations were used in an appropriate conservation form 
to insure jump conditions corresponding to the isentropic limit of the 
Rahkine-Hugoniot conditions. As a consequence of the dissipation generated 
by the truncation errors, the shock discontinuities were endowed with a profile 
(to be sure a fictitious one) , so that the flow could then be treated as 
entirely continuous. Of course, a fine mesh must be employed about the shocks 
to obtain an adequate spacial resolution of the shocks and to min imi ze the 
extent of the inadmissible region within the shock profile. Embedded shocks 
were thus evolved automatically in the analogue, a liability (truncation 
error) here being parlayed to resolve the very difficult task of treating the 
a priori unknown embedded shocks - a rare event indeed especially in transonic 
aerodynamics. A drawback of these early efforts was the resulting complexity, 
in particular in the fulfillment of the boundary condition at the airfoil, 
which led to a large expenditure of computer time. The emphasis at this 
early stage was concentrated in obtaining meaningful results, and this pre- 
occupation compounded by the conservatism engendered by the lack of both 
experience and resources resulted in a procedure, which, in light of progress 
to date, must now be considered obsolete. Thus the goal of the present effort 
is a thorough renovation of the above unsteady procedure to reduce the com- 
puter costs without significantly compromising the accuracy and the gen- 
erality of flows that can be treated. At the same time the resulting numeri- 
cal procedure must be suitable for use in an overall procedure to incorporate 
viscous effects, and for ultimate extension to steady three dimensional 
flows. The retainment of the unsteady format here to treat steady flows 
was motivated to be sure by our past experience, but mare so by the firm 
foundation on which the unsteady analogue rests, namely in the treatment of 
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shocks via the concept of weak solutions developed by Lax (Ref. 4 ). a real 
unsteady approach would, of course, also fill a dire need for a procedure to 
treat supercritical unsteady flows. 

Subsequent to the introduction of the above unsteady analogue there has been 
phenomenal progress in the development of steady finite difference procedures 
using relaxation techniques that has reduced the problem of large computer 
times. The key to this progress was the introduction by Murman and Cole 
(Ref. 5) of a finite difference scheme that was locally tailored to be com- 
patible with the local domain of dependence, and the use of a line relaxation 
in a direction natural to near-sonic flows. In their efforts the transonic 
small perturbation equation of Guderley and Von Karman was used together with 
planar boundary conditions. The basis of subsequent more exact procedures of 
Lomax and Steger (Ref. 6) and Garabedian and Korn (Ref. 7) have been centered 
on the above innovations of Murman and Cole. (These procedures are reviewed 
in Ref. 8). Despite the ostensible success of the above steady analogues, there 
are some apparent shortcomings in the treatment of the shocks that are diffi- 
cult to assess fully due to the absence of a guiding background comparable to 
Lax's contributions for the hyperbolic unsteady case. 

In Section 2 to follow we shall first as a review formulate the problem 
starting with the exact unsteady non- i sent ropic case, followed by the 
Guderley- Von Karman unsteady small disturbance case, and the small disturbance 
and exact cases with fictitious time derivatives. The intent here, in order, 
is to introduce the use of the concept of the weak solution by Lax to justify 
heuristically the' use of the unsteady analogue to treat flows with shock 
discontinuities, to review the consequences of the small disturbance limiting 
process to clarify the compatibility of shock pressure drag with the isen- 
tropic state of the flow, and finally to introduce the use of fictitious 
time derivatives to advantage when the sole interest is in the steady flows. 

In Section 3 factors in the formulation of the problem and in its numerical 
analogue which affect computer costs, will be discussed briefly, while in 
Section 4 some results using the simple equations and boundary conditions 
will be described for the severe case of a highly blunt-nosed, aft-cambered 
profile representative of the supercritical profiles now receiving extensive 
attention. In Section 5 we shall next describe the procedure proposed to 
meet the objectives listed earlier, and illustrate its use for a case with 
a "peaky” supercritical airfoil. In the concluding remarks we shall review 
and assess the treatment of the shock discontinuities by both the steady 
and unsteady finite difference analogues, (in this connection the role of 
Zierep's local analysis about the "foot of the shock" will be related); 
finally followed by a discussion stressing the need to incorporate the 
viscous effects. 
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2. Formulation of the Airfoil Problem 


We are concerned with flows over a general class of profiles at high subsonic 
Mach number ^ in which one or more embedded shock discontinuities may arise. 
The flow thus is to be considered inviscid everywhere except at the shock 
discontinuities where the flow particles are assumed to suffer an abrupt in- 
crease in entropy as they pass through the' shock. The flow is assumed to 
be non-heat conducting everywhere including the shock so that any heat ad- 
ditions by dissipation within a fluid particle is confined to that particle; 
that is, the flow will be everywhere adiabatic. The flow will be additionally 
isentropic (and therefore irrotational) at all points except in the subdomain 
downstream of the shock discontinuities, where the flow will be non-isentropic 
and therefore rotational when the shock is of variable strength or curved in. 
shape. Downstream of the shock the entropy will be constant along any given 
streamline, but will in- general differ from streamline to streamline. 

Under these conditions the basic flow equations in cartesian coordinates are 
given by 


9W + 3F + dG 
dt dx dy 


( 1 ) 


where the vectors W, F and G are given by 


pv 

r l / 2 2. 

^p[~( u + v ) + e] J 


pu 


pv 

2 



pu + p 

G =, 

puv 

puv 


2 

pv + p 

12 2 

l P- u [~( u + v ) + h] ; 


12 2 

l pv[-(u + V ) + h ]J 


4 



When the equations are given in the divergenceless form of Eq. (l) , it is 
said to be in conservation form. In the above the usual symbols have been 
employed with the less familiar symbols h and e noting respectively the spe- 
cific enthalpy and specific internal energy. For ideal diatomic gases the 
specific enthalpy is defined and given by 


h = e 



£ 

P 


while the specific internal energy e is given by 


e 



£ 

P 


( 2 ) 


(3) 


Here c . is the specific heat at constant volume, y is the ratio of specific 
heats equal to 1.4 for air, and T is the absolute temperature. With the 
value of h and e from Eqs. (2) and (3) inserted in (1), the resulting equations 
form a fully determined set of four equations for the four unknowns p, p, u, 
and v. A further useful thermodynamic relation, the caloric equation of 
state, which is in fact implied in Eqs. (2) and (3) is given by 

f - S 
— — 

v 



where R = c_ - c is the gas constant, S is the entropy, and S is a constant 
P v o 

reference value. 


Across the shock discontinuities we must additionally fulfill the shock jump 
conditions representing the conservation of mass, momentum, and energy. 

Thus in general for a moving shock with the component of the shock velocity 
normal to the shock, we have 


»i (n i 


N ) = p (N - N ) 
S h 2' 2 S 


(mas s ) 


P 1 V\ " V + P 1 = p 2 W 2^ W 2 " V + P 2 ( momentum ) 

T = T (momentum) 


(5) 


1 2 2 

T p (N - N ) (N + T ) + p (N - N ) e + N p (energy) 

2 K l v 1 S K 1 1 *1 1 S I 11 

1 2 2 

= - p.(N- - Nj (N + T ) + p (N - Nj e + N p . 

2 P 2 V 2 S ' 2 2 P 2 V 2 S 2 2 *2 
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Here N and T are respectively the velocity components normal and tangent to 
the shock front, and the subscripts 1 and 2 denote conditions upstream and 
downstream of the shock. (See Sketch A ) . 



Sketch A. The Oblique Shock Nomenclature 

In addition to the jump conditions (5) we must additionally impose the 
condition 

s 2 > s x (6) 

where the entropy S is defined in (4) . 

The jump conditions (5) can be expressed in terms of the cartesian velocity 
components by using the following correspondence: 

N. 

— = sin (a - 0.) i = 1, 2 

W. i 

l 

T. 

— = cos (a - 9 ) 

W. v i 

l 
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where 


TT 2 2 2 

W. = u. + v 
x i i 


tan 0 . = — 
x u 


and a is slope of the shock relative to the freestream direction. We note 
here for future use the well-known result that the change of entropy 
across the shock is third order in terms of the shock strength defined, 


- p^ across the shock, 
non-isentropic flow an 


for example, as the difference of the pressures p^ 

Finally we recall that for a steady adiabatic, bul 
integral of the flow equations, the strong Bernoulli law, exists valid along 
a streamline given by 


.12 12 
h + - q = - q 
2 2 H L 


(7) 


2 2 2 

where q = u + v , and q is the limit speed obtained by expanding the gas 
isentropically to a vacuum. The value of the limit speed q^ will in this 
case be constant along a given streamline but will vary from streamline to 
streamline in accordance to the equation 


12 _ — 
grad “ q = T grad S + q x curl q 
2 L 


where q is the velocity vector. Farther with the use of Eq. (2), Eq. (7) 
can be rewritten as 

P - e<i L 2 - <7> W 

We shall use this expression later to eliminate p for an isentropic case 

where q is a constant. 

E 

Consider next the formulation of the boundary and initial conditions to 
supplement the differential equations (l) for the singly connected domain 
G as enclosed by the boundary C of Sketch B. 

Here the outer boundary ABC is taken sufficiently. far from the airfoil that 
the freestream conditions can be prescribed on this contour for all times. 
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At the surface of the body y = y^(x) f° r t ^ 0 we require the flow tangency 
requirement 


v(x,‘ y ) 

£_ _ i 

u(x, y g ) y B 


(x) . 


(9) 


(Note that y (x) must be a double valued function of x, th£ two branches 
corresponding to the upper and lower surfaces 0 At the airfoil trailing 
edge we must fulfill the Kutta condition requiring the equality of pressures 
at K and D for t > 0. For an unsteady case where the lift is changing with 
time, the rear stagnation streamline will be a vortex sheet (contact surface), 
so that for t > 0 we must fulfill the jump conditions that at corresponding 
points along KA and DC the pressure and the flow direction must be equal. 

As the steady state is reached, of course, the vortex sheet disappears; and 
the flow becomes continuous across the rear stagnation streamline. At 
corresponding points on either side of the shock, HI and IJ, one must fulfill 
the shock jump conditions given earlier in Eqs. (5) and (6). Finally at a 
given time t = t < 0 we must prescribe the initial values of the dependent 
variables at all points of G. This initial flow need not fulfill any of the 
above boundary conditions, and an uniform flow at the freestream conditions 
or an available neighboring flow is generally used. 


The flow in domain G is entirely continuous, so that the flow equations may 
be used in any appropriate form; in particular, they need not be used in 
the conservation form given in Eq. (l). The difficulty of the above choice of 
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the domain G is that the location of the shocks is not known beforehand, 
or if it is known at the initial time, a procedure must then be incorporated 
to determine the shock locations at all later times, clearly a difficult 
task. A procedure based upon this approach was in fact devised ingeneously 
by Grossmann and Moretti (Ref. 9 ) who treated several simple cases. In the 
present study, we shall circumvent the task of having to determine the in- 
stantaneous shock locations by redefining the inner boundary of Sketch B such 
that it lies adjacent to the airfoil surface and does not exclude the shocks. 
The flow domain defined by such a boundary will now contain flow (shock) 
discontinuities, and the dependent variables and their first partial deriva- 
tives will no longer be continuous, but will be piecewise continuous. In 
such a situation derivatives can no longer be defined at the shock discontinu 
ities making it impossible to check the fulfillment of the equations. To 
circumvent this difficulty Lax (Ref. 4 ) used a generalization of the concept 
of a solution, the so-called weak solution, which we shall briefly review 
in the next sub -section. 


2.1 Weak Solutions. (See Courant- Hilbert, Ref. 10) 

To handle solutions of first order quasi-linear hyperbolic equations which 
are piecewise continuous, and possess piecewise continuous first derivatives, 
the concept of the weak solution was introduced by Lax (Ref. 4 ) . Thus con- 
sider the simple case of Eq. (l) in two dimensions 


W + F =0 (10) 

"t X 

Instead of requiring this equation to be fulfilled, a "-weaker" requirement 
on the "solution" is imposed given by the equation 



+ F ) dx dt = 0 
x 


( 11 ) 


where R are subdomains of G (the flow domain) , and Q are arbitrary "test 
vectors" having a certain required smoothness and vanishing outside of R. 
First for subdomains R not containing a shock discontinuity Eq. (ll) reduces 
to (10) , and the usual solution, the genuine solution, is obtained. For 
subdomains containing a shock Eq. (ll) will be rewritten as 


//C O t + f x 3 <ax dt = Jf [<fcw) t + (cf) x ] dx dt 

R X R 


" f f (WC t + FC x ) dx dt = 0 


( 12 ) 


R 
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The first integral on the right side may be reduced to a contour integral 
using Green’s formula with the result 


//COT + OT ] dx dt = <j> £[Wdx - Fdt] 


where s is the contour enclosing R. If we next use instead a contour equiva- 
lent to s shown in Sketch Q, and the fact that C. = 0 along those segments 
of the alternate contour not "wetting" the shock, we can rewrite ( 13 ) as 


/ /[OT t + OT X J d* dt = / C,[W] 


shock 


— - [f] — 

a . - as 

shock shock 


where ds is a line element along the shock, [w] 3 ^ - W g and [F] = F^ - F g 
(subscripts 1 and 2 denoting conditions upstream and downstream of the shock), 

dx dt 

and — and — are the direction cosines of the shock, 
ds ds 

shock shock 

Eq. (12) using Eq. (l4) now becomes 


/ f(W£ + F£ ) dx dt - f Qi 


ds = 0 


shock 


shock 


shock 



ORIGINAL CONTOUR FOR R: ABGHA 

REVISED CONTOUR: ABCDEFGHA 

Sketch C. Revised Contours for the Weak Solution 
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We now define a weak solution W (or F) of (10) if it is piecewise continuous 
with piecewise continuous first derivatives, and if it fulfills the equation 

//(we t + F £^) dx dt = 0 ( 16 ) 


for all admissible "test vectors" £ in all subdomains R within the flow 
domain G. From (15) with (16) fulfilled, it follows that weak solutions 
will fulfill the condition 


/«• 

_ dx 

rW1 
' [F] 


shock 


shock 

shock 


(17) 


Since £ is an arbitrary "test vector", for Eq. (17) to be fulfilled we must 
require 


[W] 


dx 

ds 


rnl dt 

' [F] a! 

shock 


= 0 
shock 


( 18 ) 


which is, without surprise, the unsteady one dimensional shock jump conditions. 
Thus a weak solution as defined above will fulfill the correct jump conditions 
across embedded shocks. The definition of weak solutions as given above can 
be extended to higher space dimensions. (See Ref. 10). Weak solutions as de- 
fined above have some surprising properties, and these are described by Lax in 
Ref. 4. Two of these are of particular relevance to us. The first is the fact 
that the shock jump conditions will depend upon the specific conservation form 
of the basic flow equations used as would be expected by inspection of Eq. (l8). 
Thus it is essential that the basic flow equations be used only in the form 
given in Eq. (l) for subdomains containing a shock. The second property is 
that for an initial value problem, the initial data alone do not determine an 
unique solution. A supplementary condition must be additionally imposed, but 
the nature of the additional condition is known unfortunately only for a few 
simple initial value problems. In a number of these simple cases with a 
first order quasi-linear hyperbolic equation where the exact analytical 
solution is known, it has been shown analytically using (l6) that an unique 
weak solution was obtained as the limiting solution (viscosity -*0) of a cor- 
responding viscous problem where an artificial viscous term was added to the 
original equation. Further Lax has shown that for the same problem the same 
weak solution was essentially obtained as a limiting solution for small mesh 
using an appropriate finite difference analogue. This is a fortunate result, 
though not altogether surprising, since an artificial viscosity is also 
implied in the first order diffusing difference scheme used by Lax. 
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For the hyperbolic combined initial-boundary value problem of higher dimension 
of interest to us, we have no direct prior analysis to guide us as in the 
simpler cases; we must extrapolate the conclusions from the simpler cases 
studied by Lax and others, and indeed we must be grateful that such infor- 
mation is available, since it does offer a reasonably sound point of departure. 
Thus we may conclude that a reasonable procedure for the airfoil problem is to 
pose it as a hyperbolic initial-boundary value problem using the equations in 
an appropriate conservation form and then solving it by an appropriate finite 
difference analogue. In such a procedure the correct shock jump conditions 
can be expected to be fulfilled across the embedded shocks. It is essential 
here that in subdomains containing a shock discontinuity, the proper con- 
servation form of the equations be used. Lax (Ref. 4) and Ref- 11 have 
shown that for the one dimensional unsteady shock wave problems when the 
equations are used in the non-conservation form, either the wrong shock speed 
results when the shock is moving, or the required jump conditions are not 
fulfilled when the shock is stationary. In passing one should finally note 
that a comparable background and point of departure does not exist in the 
steady finite difference analogues. 

2.2 Guderley-Von Karman Formulation of the Problem (See Guderley, Ref. 12) 


We need not dwell at length on the transonic small disturbance theory except 
to recall that it is the result of a very definite limiting process t -*0 
carried out in the basic equations (l) using the transformation (hypothesis) 


u(x,y,t) = a*[l + T u(g,Tl,?')] 
v(x,y,t) = a* t 3//2 v(g,Tl,t) 
P(x,y,t) = p*[l + t 3 ^ 2 ■?(§,T1,'£)] 
p(x,y,t) = p*[l + t 3 ^ 2 ?(S>T),t)] 


where 



„ 1/2 
T\ = t 

X 

o 


t = U)t 


(19) 


For a limit of to representing reduced frequencies of order unity, the result- 
ing limiting equations (lowest order in t) are found to be 


uL = -(y +l)uu+v ; 0 = u - v 

t 5 11 T) ? 


( 20 ) 
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Applying (19) to the airfoil + boundary condition, one obtains the standard 
planar condition v(§ ±0) = h (§) where h - (g) are the slopes on the upper and 
lower sides of the profile for the stretched coordinates. There is mu ch 
implied in Eq. (19)5 and in particular the various powers of t have of course 
been selected here, so that the resulting lowest order equations and the 
coordinate stretching would be physically viable. The parameter t (propor- 
tional to the airfoil thickness ratio to the two-thirds power) is a measure 
of the proximity of the flow velocities to the sonic speed, so that the limit- 
ing process t -*0 implies values of Mco close to one and thin airfoils at small 
incidences. It is particularly interesting that in this limit the flow is 
irrotational and hence isentropic. This should not be surprising since the 

p - p* 

pressure coefficient C s is first order in x 5 while the entropy 

P q* 

3 

change across the shock is third order in the shock strength, or 0(t ). 


The order of magnitude assessment of the various terms in the flow equations 
solely in terms of powers of r can be misleading since the actual magnitude 
of a term in fact depends additionally on the algebraic magnitude of 'the 
term itself. Thus in an unfavorable situation it may require such a small 
value of x for the above small disturbance theory to be valid, that the 
resulting flow would be of no practical value. 


Lastly, it should be noted that the unsteady Guderley-Von Karman equations 
as given in Eq. ( 20 ) are of mixed parabolic-hyperbolic type, so that Lax’s 
concept of weak solutions does not apply directly to these equations. If 
we are interested only in the steady state, fictitious fully hyperbolic 
equations can be obtained by introducing appropriate unphysical time deriva- 
tives, and we shall examine this in the next sub-section. 


2.3 Fictitious Unsteady Equations to Treat Steady Flows 

In preparation for Sections 4-5 in the present sub-section we shall simply 
note several sets of equations all under isentropic conditions, having 
fictitious time derivatives which result in fully hyperbolic systems when 
added to either the Guderley-Von Karman or the exact steady terms. . Thus for 
the small disturbance case we have 


U~ = - (y + l) (~ U ) + 

t \ 2 \ 


(v) 



v 

§ 


( 21 ) 
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or with the exact steady terms (in non- conservation form; 


2 2 2 2 
u = (a - u ) u - 2uv ii r +(a - v ) v 


v = u - v 
t y x 


where the speed of sound is given by 


2 1 2 1 2 2 
a = g(y + 1) a* - -(y - l) (u + V ) 


Eq. (22a) can be also rewritten in a conservation form 


u + = (p u ) v + (pv) 

u x y 


V = - V + u 
t X y 


where p = p(u, v) from the isentropic relation. 


(22a) 


(22b) 


It is to be noted that these equations are not in the proper conservation form 
that yields the correct shock jump conditions, and its use should be restricted 
to subdomains where the flow is continuous. 


Lastly we note for subsequent use the exact unsteady Euler equations in the 
proper conservation form obtained from the Eqs. (l) for the isentropic case, 
that is, 


St 


where 


r ' 


f ' 

p 


pu 

pu 

S_ 

+ dx < 

2 

pu + p 

pv 


puv 

. J 


k J 



r i 

pv 

, P uv . 
2 

pv + p 


= 0 


E_ 

P* 


/ JL.V = .Y-T \ £_ 

l p* J 2y p* 



( 23 ) 


Shock jump conditions for the asymptotic steady state for the above hyperbolic 
equations can be obtained using Eq. (l8) extended to the unsteady planar flow 
case. Thus for the small disturbance equations we obtain the jump conditions 
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[- i( Y + 1) U 2 ] f 

L 2 Y dx 


- [£) = 0 




i 


dx 


shock 

+ [u] = 0 


shock 


Here again the brackets [ ] denote a jump of the bracketed quantity across the 

shock. Thus for a normal shock adjacent to the airfoil surface, where [v] = 0 

in accordance to the planar airfoil condition, we have [u ] = - u 0 ^ = 0, 

or u = - u > since u f u . 

1 2 1 r 2 


Such a condition conforms to the transonic shock polar only if v ~ 0 , that is 
for shock waves essentially normal to the free stream. In the planar case if 
the basic small disturbance restrictions are not violated, shock waves that 
depart from a normality with the free stream cannot arise, so that the shock 
captured by the small disturbance weak solution will automatically conform to 
the consistent shock polar. The same situation would prevail also for the 
three dimensional case described by the usual small disturbance equations, 
where again shock jump conditions implied by the weak solution will not conform 
to the appropriate shock polar if the shocks are not essentially normal to the 
unperturbed basic flow direction. Thus if an excessive wing sweep causes an 
excessively swept shock (sweep of 25° is excessive), then the small disturbance 
jump conditions would not be fulfilled for such shocks. (Excessive sweep effects 
have of course never been permitted by the small disturbance hypothesis.) 


Lastly we note that the shock jump conditions for the exact equations (23) in 
the asymptotic steady state are given by 
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which after some algebraic manipulations can be shown to agree with the 
Rankine-Hugoniot conditions given earlier in Eq. ( 5 ) with the shock speed 
set equal to zero. 


2.4 The Paradox of Isentropic Transonic Wave Drag 

It is universally accepted that the onset of transonic wave drag over airfoils 
in an inviscid flow for < 1 is the consequence of streamwise momentum loss 
in the far field due to entropy additions by shock waves. The question then 
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naturally arises as to how a drag can arise if an isentropic flow is assumed. 

To answer this let us consider the determination of the drag by the vanishing 
of the flux of the total streamwise momentum (including the pressure term) 
through the control surface of Sketch B. Consider the case that the flow is 
calculated with the exact non- isentropic equations with the shocks treated either 
as discontinuities fulfilling the Rankine-Hugoniot jump conditions 3 or evolved 
by the continuous shock-capturing weak solution to the appropriate conservation 
equations . Vanishing of the streamwise momentum flux through the contour C 
(Sketch B) then yields the equivalence of the airfoil drag with the deficiency 
of the streamwise momentum flux through the far downstream portion of the con- 
tour, since the momentum is conserved across the shock. Since this far field 
momentum deficiency can be directly attributed to the increase in entropy, in 
the exact non-isentropic case the drag is indeed the consequence of the shock 
entropy additions. If we now calculate the flow in the same manner but using 
instead the exact isentropic equations and evaluate the total momentum along 
the control surface, we find that the momentum flux deficiency along the far 
downstream boundary vanishes identically, and the drag is now equal to the 
change of the streamwise momentum flux across the shock wave. If the equations 
are used in the proper conservation- form about the shock (as is done in the 
present report) , then the momentum is conserved across the shock,- and the drag 
must then be identically zero. Care must be therefore exercised in the choice 
of the proper form of the conservation equations, since if the equations are 
used in the improper conservation form or in the non-conservation form such 
that the conservation of momentum is not preserved across the shock, then an 
erroneous pressure drag would result precisely equal to the erroneous decre- 
ment in the momentum conservation across the shock. In a numerical calcula- 
tion truncation errors seldom permit the drag to vanish exactly; the finite - 
ness of the drag would then be a measure of the accuracy of the calculations. 

If, however, the Von Karman-Guderley transonic small perturbation equations 
are used, again in such a way that the shock jump conditions in a consistently 
simplified form are fulfilled, one finds that an exact integration of the 
streamwise component of the surface pressures will not integrate to a zero 
drag. This apparent paradox results since the conservation of momentum across 
the shock is fulfilled only to within the small disturbance approximation, 
and the drag that is obtained would just match the decrement in momentum 
conservation due to the inadequacy of the simplified theory. (Provided of 
course the nose region is calculated properly, which is seldom accomplished.) 

Thus in a properly formulated isentropic theory the consideration of drag, in 
particular, drag divergence, is beyond the scope of the theory, while on the 
other hand, the determination of the pressure distribution, and hence lift 
and moment, is well within the prediction scope of the isentropic theory so 
long as the resulting shock strengths are no stronger than those experienced 
in experiments. 
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The above remarks apply of course only for < 1. For = 1 and > 1 
shock pressure drag can arise in an isentropic flow by the loss of streamwise 
momentum by radiation through "supersonic windows" that open up in the lateral 
far field. 


3- FACTORS AFFECTING COMPUTER COSTS 

Consideration of the computer costs in a finite difference analogue must 
already commence with the manner of formulation of the initial-boundary value 
problem on which the analogue is based. The equations should of course be 
used in the simplest form possible , both in terms of the complexity of the 
terms and their numbers as well as in the number of equations itself, weighed 
of course with the desired range of applicability of the resulting equations. 
Thus the simplest meaningful set of equations in the unsteady approach is • 
the Guderley-Von Karman equations with fictitious time derivatives, Eqs. (2l). 
These equations are in the proper conservation form conforming to the defi- 
nition of Lax's weak solution. Here the assumption of small disturbances will 
impose limits on the class of profiles that can be treated. In the steady 
case with the small disturbance equations one has the choice of either having 
two first order equations in u and v, or one second order equation in the 
perturbation potential. Here the reduction of the number of equations must 
be weighed against the increase in the order of the equation, as well as the 
increased accuracy with which the potential must be determined to obtain a 
given accuracy in its derivatives which are the physically meaningful 
quantities. Having to deal with one quantity, the potential, instead of the 
two velocities will, on the other hand, reduce the computer storage require- 
ments. Simplest equations frequently arise when they are expressed in 
cartesian coordinates. Since the latter coordinate system invariably leads 
to a near-intolerable incompatibility with the boundaries on which conditions 
are prescribed, some compromise here must be sought, though generally the 
simplicity of the equations must be sacrificed at least locally. 

Consider next the boundary conditions, in particular that at the airfoil 
surface which is probably the most crucial. With respect to the condition 
itself, it is invariably most preferable to prescribe an explicit non-homo- 
geneous condition in terms of only one of the dependent^variables; that is, 
it is more preferable to have a condition such as v = h~(x) at the airfoil 
surface, than v - h u, where h is a known function. The "tieing down" 
of the value of one of the unknowns at the boundary here should in a numerical 
procedure which is marginally stable lead to greater stability than if the 
homogeneous condition were used. Where only the more, exact homogeneous 
airfoil condition, suffices, it will be more effective to set initially 
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u = u_ (freestream value) to obtain initially an explicit non-homogeneous 
condition of the type v = h~(x), and then upgrade u by an iterative procedure. 
In imposing the condition it is most convenient to be able to impose it, 
not on the exact boundary, but on a neighboring boundary which is more com- 
patible with the coordinate system, as in the case of the classical planar 
conditions. Where such planar conditions are inadequate, as we shall soon 
see is the case in the nose region of a blunt-nosed body, it may be sufficient 
to fulfill the exact conditions on a simplified neighboring non-planar boundary 
conforming exactly to the airfoil only where it is necessary. We shall 
illustrate one example of this in a later section. Another example where 
such a simplification might be used is in the procedures where the airfoil 
is mapped onto a circle. Here a simpler mapping might be considered, where 
only the nose region for example is mapped exactly onto the circle while the 
remainder of the profile .is mapped to the neighborhood of the circle. The 
exact boundary condition would, of course, be prescribed on the circle in a 
"quasi -planar" way. 

Let us turn next to the finite difference analogue for the unsteady case 
assuming here that we start with the simplest formulation of the problem 
suitably compromised to be consistent with our requirements. There are three 
factors that now determine the overall computer time. The first factor con- 
cerns the amount of calculation required for an elemental step where the 
solution at a given point is advanced to the next time step. Here an im- 
portant consideration is the choice of the difference scheme. There are of 
course a multitude of difference schemes which have been suggested, but it 
will suffice here to discuss explicit schemes only, and among these only 
Lax's first order diffusing scheme (Ref. 4), and the two step Lax-Wendroff 
second order scheme 'described by Richtmyer (Ref. 13). There have been 
suggested other more elaborate schemes of second order accuracy which appear 
to be more effective under special conditions, but in general the added 
complexity beyond the straight-forward Lax-Wendroff scheme does not appear 
to be warranted. Quite obviously the second order Lax-Wendroff scheme for a 
given mesh spacing will yield a smaller truncation error in the marching 
process, and. hence a more exact solution than a first order scheme, but of 
course at the expense of requiring more arithmetic steps. In the problem 
in the (x, y, t) system, in the second order scheme five first order schemes 
are in fact required. The mesh size required in the x, y, or t directions to 
obtain a given accuracy depends on the gradient of the dependent variables 
which arise in. these directions. Assume that in a given case that the required 
mesh dimensions Ax and Ay are equal. The largest step At that is permitted 
for linear stability is dictated by the Courant, Friedrichs, and Lewy con- 
dition in the form, At < const. Ax, a statement requiring the point at the 
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new time to fall within the region of influence of the initial data. Thus 
in a hypothetical case assume that with the first order scheme half the mesh 
size Ax used for the second order scheme is required to obtain the same 
accuracy. This then means that 2^ elemental calculations must be carried out 
with the first order difference scheme in contrast to one elemental step 
involving five first order differences for the second order scheme, a reduction 
of computer time by a factor 5/8 for the latter scheme. Unfortunately in a 
given unsteady planar flow problem the situation is not so clear cut, with 
the precise merits of the first and second order schemes for a sufficiently small 
spacial mesh being determined largely by the time rate of change of the develop- 
ing flow. Quite clearly if the flow is only changing slowly with time, the 
first order scheme may be entirely adequate, while on the other hand the first 
order scheme will be entirely inadequate when the flow has experienced a large 
unsteadiness prior to its approach to the steady state. In the latter situa- 
tion the first order scheme is quite adequate to handle the slowly changing 
approach to the steady state, but it has suffered a large truncation error 
during the unsteady phase from which it can never recover, resulting in an 
erroneous steady asymptotic flow. The latter behavior has been found in 
several calculated examples involving a blunt-nosed supercritical profile. 

In fact for these test cases the first order scheme was found to be entirely 
inadequate, and only the use of the second order scheme led to viable results. 
Quite obviously an attractive procedure suggested here is the use of the 
second order scheme to treat the more unsteady phase of the time history, 
and then switch to the first order scheme in the slow approach phase to- the 
steady state. 

In passing it may be noted that the use of the two step second order Lax- 
Wendroff difference scheme in the unsteady planar equations in conservation 
form leads to a "telescoping" feature leading to the automatic fulfillment of 
the integrated conservation laws. (See Richtmyer, Ref. 13)* 

The second factor influencing the computer time is the total number of mesh 
points which are required ultimately in the x, y plane to describe the final 
steady flow. The mesh spacing here is simply dictated by the flow gradients. 
Quite obviously economy would suggest the use of a variable mesh spacing to 
fit the final velocity gradients, but this cannot be overdone since a waste- 
ful overlap of the different sub-domains with a given mesh must be provided 
to assure continuity of the flow. In the early efforts it was thought that 
the initial phases of the developing flow could be treated by a coarse 
spacial mesh to establish the intermediate and far field flows rapidly, and 
only after this introduce a refinement of the mesh about the airfoil to 
recover a viable near field flow. This simply did not work out, at least for 
blunt nosed supercritical profiles, since the strength of the "singularities" 
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to be communicated to the intermediate and far fields were never adequately- 
developed in the near field, due to the coarse mesh, in particular, the 
contribution of the nose flow. Another means to decrease the number of mesh 
points is simply to redefine the dependent variables such that the sought 
quantities at the new time are differences from the initial values. Such 
differences will have significantly decreased gradients, especially in the 
long duration approach to the steady state, so that the spacial mesh size 
can be accordingly increased. Since as we have seen earlier, the computing 
time varies inversely as the cube of the mesh dimension, such a simple 
change of dependent variable should have a dramatic impact especially in the 
vicinity of a highly blunt nose where generally an early attainment of steady 
conditions results and where a disproportionate number of mesh points is 
required. In this change of variables, most probably a linearization can be 
justified except perhaps about the shock if the shock location is changing. 

In a version suggested in Ref. 14 to be .studied shortly, the non-linearity 
was retained. Here it is important to reexamine the impact on the shock 
jump conditions. 

Finally the third factor affecting computer costs is the rapidity with which 
the desired steady state is obtained. The last procedure described above 
will clearly have a bearing here also. Another important facet here is the 
choice of the initial data itself. Quite clearly the proximity of the initial 
flow to the final flow is of significance, but the route between the two is 
not always direct as was found in several calculated examples , where a slowly 
oscillating apparently near-undamped shock sequence was obtained. The simple 
cure here is simply not to prescribe an initial flow which is susceptible to 
this phenomenon (though we do not of course have a clear cut criterion for 
this) by using "adjusted" initial data which places the shock somewhat up- 
stream of the expected location to provide a fast adjustment path for the flow 
to adjust to the desired state. 

In summary the path' suggested to reduce computer time without sacrifice of 
the accuracy of the solution is to tailor not only the formulation of the 
problem but also its analogue, both in space and time, to be "just adequate." 
Thus in regions not containing a shock the equations in the simpler non- 
conservation form is suggested, while about the shock the appropriate con- 
servation forms must be used. At a blunt nose the exact boundary condition 
should be applied at the exact surface while elsewhere the condition may be 
fulfilled in a quasi-planar way on a simpler boundary in the neighborhood of 
the airfoil surface. In the analogue a variable spacial mesh spacing to fit 
the anticipated flow gradients should be used; the two step second order Lax- 
Wendroff difference scheme must be used during the initial unsteady phase 
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converting to a first order scheme with a change of the dependent variables 
without a linearization in the final time consuming approach to the steady 
state. Quite obviously the extent to which the above philosophy should be 
implemented must be weighed against the increased program logic and storage 
required as well as by the wasteful overlap and added complications to avoid 
a mismatch at boundaries where the procedures are patched. 

In the next section we shall first examine whether the Guderley-Von Karman 
equations with fictitious time derivatives and with planar boundary con- 
ditions are adequate to treat blunted nosed supercritical profiles. This is 
then followed by a more exact formulation incorporating some of the suggestions 
given above. 


4. INVESTIGATION OF THE SIMPLIFIED PROBLEM 

In the present section we shall first describe the finite difference results 
using the Guderley-Von Karman small disturbance equations as given in Eq. (2l) 
with a planar boundary condition at the airfoil. The test case will be the 
flow at M^ = 0.80 over the supercritical profile of 11$ thickness ratio at 
a = 0 calculated earlier in Ref. 3 with the exact conservation equations (23). 
This case was selected because 

a. It is fairly representative of the current problems on which one 
might want to employ the numerical analysis. 

b. The moderately large thickness (11$), heavy cambering, considerable 
nose blunting (2.2$) and fairly high local Mach numbers would provide 
severe tests of the ability of the simplified equations and boundary 
conditions to describe adequately the flow. 

c. A solution obtained by the numerical method described in Ref. 3 was 
available for comparison. 

In the calculations a modification of Richtmyer's two step version of the 
second order Lax-Wendroff difference scheme has been used with an inclusion 
of some diffusive damping. (See Ref. 2). As is well known this explicit 
difference scheme is stable only if the marching increment At is kept below 
some limiting value dependent on local conditions. For simplicity the allow- 
able At has been determined for a representative range of u and v to be en- 
countered in typical airfoil problems. The complex amplification of small 
disturbances of various wave lengths and orientations was determined as a 
function of the marching increment. At, and limiting values which would keep 
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the absolute value of the amplification below unity were chosen. These 
limits were then examined as functions of u and v, and an overall limit, 

At < 0.15 Ax was chosen. This limit is highly conservative around the stag- 
nation point and is closest to the actual stability limit for sonic flow at 
an inclination near 30° to the x direction. 

The preceding remarks on stability and damping apply to the centered differ- 
ence schemes used at ordinary field points. At the airfoil surface (or, 
rather, the chordline slit) where boundary conditions are to be imposed, 
compatible non-centered difference schemes were devised. That is to say, 
special schemes were constructed to have the same nominal order of accuracy 
as the ordinary schemes. These were also screened by checking amplification 
of small disturbances of various wave lengths and orientation. . This screening 
was considered necessary but insufficient to prove stability inasmuch as the 
presence of the boundary can shift disturbance energy into disturbances of 
different wave length and orientation; this would be incompatible with the 
fundamental assumptions in the linearized stability analysis which is being 
employed. The boundary 'schemes are of the Lax-Wendroff one-step class with 
special matrices for the stabilization terms. The orientation of the boundary 
with respect to the mesh influences the magnitudes of the elements of the 
stabilization matrices. 

In the calculations a cartesian mesh system was employed about the airfoil 
overlapping with a polar mesh system in the far field. Approximately 3500 
mesh points were typically used distributed as shown in Table 1. Here the 
airfoil chord is of unit length with the nose located at the origin. 


In the present problem the flow about the leading edge is singular. To avoid 
this region in the calculation, a similarity solution for a parabolic nose at 
zero incidence is used in a 0.0125 chord square domain centered about the 
slit terminus. This similarity solution is defined in terms of the perturba- 
tion potential by the hypothesis 

9 = hl u/7 f(c) 

where 

- 6/7 

c = (y + 1) i h 1 > and 

(36£ 2 - 49F') F" + 30£F' - 12 F = 0 . 
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Table 1 . Typical Mesh Arrangement 


Region 

Ax = Ay 

Range 

Inner Nose Region 

0.00625 

-.075 < x < .075 

-.10 < y < .10 

Outer Nose Region 

0.0125 

-.15 < X < .15 
-.20 < y < .20 

Airfoil Region 

0.05 

-.40 < y < .60 

Pine Outer Region 

0.10 

-.70 < x < 1.70 
-.80 < y < 1.20 

Coarse Outer Region 

0.20 

-1.10 < x < 2.10 

Outer Polar Region 

Ar > 0.2 
A6 = 9 ° 

1.4 < r < <a 
-9° < e < 360° 


In Fig. 1 we now show the resulting surface pressure distribution in terms 
of the pressure coefficient C given in accordance to small disturbance 
theory as ^ 



Here for comparison is shown the results from the more exact calculations 
.from Ref. 3. The comparison shows a drastic discrepancy in the nose region, 
both on the upper and lower surfaces, with the simplified solution exhibiting 
a horrendous nose over -expans ion. There is further a discrepancy in the shock 
location and in the over pressures on the lower aft surface. In the plateau 
region, on both the upper and lower surfaces, there is good agreement. It 
is thus amply clear that the small disturbance simplifications are inadequate 
at least for the test example. 

The cause of the above discrepancies can be attributed to the combination of 
the inadequate equations, and an inadequate airfoil boundary condition imposed 
at an inappropriate boundary. To separate these effects we have next calcula- 
ted for the same test flow the following four additional cases characterizing 
each case in the following by noting the equations and the airfoil boundary 
condition used. 
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Figure 1. Results of the Guderley-Von Karman 
Small Distribution Problem 
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Figure 2(a). Mesh Arrangement for Calculations Based on Blunted Slab 
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■Figure 2(b). Mesh Arrangement for Calculations Based on Blunted Slab 
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I. 

Equations : 
B.C. : 

Exact non-conservation form - Eq. (22) . 
Planar condition. 

II. 

Equations : 
B.C. : 

Small disturbance. 

Exact condition [ v = (l + u) tan 9] fulfilled on a 
slab described earlier in a first order planar fashion. 

III. 

Equations : 
B . C . : 

Exact in non-conservation form. 
Exact (Same as II). 

IV. 

Equations: 
B.C. : 

Exact in non-conservation form. 

Exact condition fulfilled on slab in second order. 


For the cases where the airfoil condition is fulfilled on the slab, the mesh 
configuration is shown in Figure 2 and detailed in Table 2. The resulting 
surface pressure distribution for these cases are shown in Figs. 3-6. Con- 
sider first Fig. 3 (Case i) with the exact equations with planar conditions. 

It is seen here that the nose overexpansions have been somewhat subdued on 
the upper surface by the more exact equations, but have been "over-subdued" 
on the lower surface to the extent that the plateau pressures have been 
additionally distorted. Planar conditions are clearly inadequate in the nose 
region. The previous discrepancy in the overpressures on the lower aft surface 
has been remedied by the more exact equations, so that planar conditions 
appear to be adequate except for the nose region. The shock wave is still 
inadequately treated, but this is most probably the consequence of not using 
the equations in conservation form. 

The results of Fig. 4 (Case II) with the approximate equations and the exact 
boundary condition illustrate the importance of having the proper condition 
at the nose. The discrepancy due to the inadequate equations at the nose 
now appear to be reduced to the order of that along the aft portions of the 
airfoil. Fig. 5 (Case III) shows that the use of the exact equations in non- 
conservation form and the exact slab boundary conditions gives reasonable 
agreement with the reference calculations except at the nose where a slightly 
increased loading arises, and at the shock. In Case IV an improvement of 
Case III is sought by simulating the fulfillment of the boundary condition 
at the actual airfoil surface by imposing a modified condition on the slab 
(second order planar condition) which assumes the velocities u and v to vary 
linearly from the actual airfoil surface to the slab. The results for Case IV 
in Fig. 6 show some Improvement due to the improved boundary condition, but 
there still remains a small discrepancy on the lower surface. 
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Table 2. Mesh Arrangement, Blunted Slab Approximations 


Region 


Range 

Cartesian 1 

Ax = AY = 0.2 

-1.1 < X 
-1.525 < 

Cartesian 2 

Ax = Ay = 0.1 

-0.7 < x 
-0.725 < 

Cartesian 3 

Ax = Ay = 0*°5 

-0.4 < x 
-0.425 < 

Cartesian 4 

Ax = Ay - 0.025 

-0.15 £ : 
-0.175 < 

Cartesian 5 

AX = Ay = 0.0125 

-0.1 < x 
-0.125 < 

Cartesian 6 

AX = Ay = 0.00625 

-0.05 < 
-0.075 < 

Curvilinear Nose 

A© = 4.5° 
0.0040 < Ar 

90° ^ e 
0.025 < 

Outer Polar 

Ar > 0.2 
— 0 

A9 = 9 

1.4 < r 
0 ~ „ 

-9 < e 


< 2.1 

y < 1.675 

< 1.7 

y < 0.875 

< 1.4 

y < 0.575 

x < 0.2 

y < 0.175 

< 0.15 

. y < 0.125 

x < 0.1125 
; y < 0.075 

o 

< 270 

r < 0.0625 

< =0 

< 360° 
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Figure 4 . Results for Case II - Small Disturbance 
Equations and Exact (First Order) Slab 
Conditions 









Thus in summary we note the following: 


b. The exact condition must be prescribed at the actual airfoil surface 
in the nose region, while it may be prescribed in a quasi-planar 
manner over the remainder of the profile. The slab condition should 
thus be adequate. 

b. The exact equations must be used, and they probably must be used in 
conservation form about the shock. (We shall check the latter in 
the next section. ) 

It must be remembered that these conclusions pertain to the more "severe" 
profiles and flow conditions, and that the Guderley-Von Karman formulation 
of the problem is most certainly adequate for a wide class of other less 
severe airfoils and flow conditions. 

5. PROCEDURE .WITH THE PROPER EMBEDDED CONSERVATION EQUATIONS 

In the unsteady procedure for further consideration, the equations with the 
exact steady terms in u and v with fictitious time derivatives as given in 
Eq. (22) will be used in the continuous regions of the flow; and the exact 
unsteady Euler equations in p, u, and v in the proper conservation form as 
given in Eq. (23) will be used in the subdomain embedding the shock. There 
is then an incompatibility of the time derivatives which would vanish as the 
steady state is attained. The exact boundary condition, as previously, will 
be fulfilled on the slab. In the subdomain embedding the shock a fine mesh 
of 1.25% chord will be employed. 

With the above system of equations several examples have been calculated. 

Two of the examples were calculated earlier with the original unsteady pro- 
cedure of Refs. 2 and 3> and they will serve as check cases for the revised 
procedure . Additionally in the present calculations reduced diffusive damping 
will be employed to capture the shocks more exactly. 

The first example calculated is the 11% aft-cambered profile at M =0.8 and 
a = 0 which was computed earlier in Ref. 3 (see also Ref. 8). The resulting 
surface pressure distribution is compared with the earlier result in Fig. 7- 
Here the agreement is reasonably good over most of the airfoil, but there 
still remains the discrepancy in the nose pressures observed earlier, and the 
shock has now been displaced essentially to the trailing edge, where an in- 
sufficient flow resolution with the mesh used precludes its proper capture. 

The obvious suspect causing these discrepancies is the application of the 
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airfoil boundary and Kutta conditions on the slab rather than on the actual 
airfoil. To check first the extent that the streamline slopes 0 differ at 
corresponding points on the slab and the airfoil, we have shown in Sketch D 
the lines of constant 0 (isoclines) from the earlier calculation about the aft 
portion of the profile where the largest differences can be expected. It is 
of interest here to note first the differences of the isocline patterns in the 
different subdomains, but only in the subsonic region aft of the shock are the 
isoclines such that a significant difference arises in the 0 distributions 
along the slab and the airfoil surface. In the present calculations we have 
prescribed the actual airfoil slopes along the slab, whereas Sketch D would 
indicate that less inclined slopes should have been prescribed. The present 
results with the slab would thus represent the flow for a profile with a some- 
what greater aft camber than desired, and this difference could readily explain, 
the greater expansion at the nose and the more downstream location of the shock. 
To verify this we have repeated the calculations using a more exact planar 
condition applied on the slab which takes into account the variation of 0 
between the slab and the airfoil. The second order planar condition used 
assumes a linear variation of 0 with y in this neighborhood where the gradients 
0 are taken as the- values obtained in the calculation at the slab. The re- 
sulting pressure distribution is shown in Fig. 7. The shock here has now been 
displaced upstream of the trailing edge and has been properly captured, but 
there now arises a small discrepancy in the plateau pressures on both the upper 
and lower surfaces indicative of insufficient aft camber. The probable cause 
of this deficient aft camber effect is the inadequacy of applying the Kutta 
condition on the slab. Examination of the earlier more exact results, where 
the boundary and Kutta conditions were imposed on the airfoil, shows that the 
pressures at the points on the slab corresponding to the trailing edge are not 
equal, but in fact a small difference arises characteristic of a case with a 
jet flap. Thus to regain the proper plateau pressures one must now permit a 
small jump in the trailing edge pressures on the slab, where the jump can be 
ascertained in the same manner as the slopes 0 were previously corrected. 

In the calculations where the first order planar slab conditions were used, as 
noted earlier the plateau pressures agreed closely with the earlier exact case, 
but in light of the above discussion this close agreement is somewhat fortuitous 
since the effects of the deficient first order planar condition apparently were 
just compensated for by the missing jet flap effect. Thus for extreme profiles 
as the present one, it is essential to apply on the slab the second order 
corrections to the slopes, while less important, a correction should also be 
considered for the trailing edge pressures. 

The second example calculated is that for the NACA 64 a 410 profile at = 0.72 
and a = 4°. The results of the present calculation are compared in Fig. 8 
with a number of previously calculated pressure distributions from several 
methods as well as with the experimental results of Stivers. Directly com- 
parable here are the present results, the result from Ref. 2, and the results 
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Sketch D. 


Isoclines for the 11 
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of Steger and Lomax (Ref. 6) . Consider first the comparison of the results 
from the unsteady procedures. The agreement of the pressure distributions 
over the lower surface is good, but there is a puzzling not insignificant 
discrepancy in the upper surface pressures ahead of the shock, most probably 
attributable to a discrepancy in the nose expansions. To check whether the 
use of the slab might account for a difference in the nose expansion, the 
fulfillment of the surface tangency condition at the actual airfoil nose 
was checked in the present results. In this highly inclined segment of the 
airfoil the resulting slopes at the actual airfoil surface were found to 
agree with the actual slopes to within less than a tenth of a degree, so 
that the use of the slab most probably must be ruled out as the source of 
the upper surface discrepancy. A comparison of the isobar distributions in 
the region about the nose however showed a significant difference. Whereas 
in the present calculation, where a more natural polar mesh system is used, 
the isobar pattern is smooth and near circular in shape, the isobars in the 
earlier calculation with a cartesian mesh system exhibited some irregularities 
(see Fig. 5, Ref. 2), and they were less reasonably shaped in general. One 
must therefore tentatively conclude that the difference of the mesh con- 
figurations about the nose and the patching of the various regions of different 
mesh sizes most probably were the causes of the difference in the nose ex- 
pansions. The results of the present calculations, despite the fact that 
they now depart from Stiver's experimental results on the upper surface 
should supersede the earlier results. (Note that the early experimental 
results of Stivers, obtained in a closed tunnel at a low Reynolds number 
without a boundary layer trip, may be afflicted with wall interference effects 
as well as more importantly, viscous effects which may cause the lowering 
of the upper surface suctions ahead of the shock through a highly probable 
leading edge separation bubble) . The use of the first order planar slab 
condition here will of course affect the upper surface pressure plateau, 
but the consequence should be relatively insignificant since the gradients 
of the streamline slopes are negligible about the aft portion of the upper 
surface. The present result further exhibits the proper jump condition, 
and there is sufficient resolution to yield a post-shock cusped pressure 
bucket. The improvement in the fulfillment of the jump conditions is also 
shown in Fig. 15 . The present results for the NACA 64A 410 profile is also 
directly comparable with the results of Steger and Lomax shown in Fig. 8 
where also the full steady Euler equations were used. The discrepancy between 
these results on both the upper and lower surface apart from the shock loca- 
tion and pressure rise, appears to be due to a difference in circulation 
brought about by the difference in the manner of imposing the Kutta condition. 

The third example is also for the flow over the NACA 64A-410 profile but at 
M® = 0.735 and a = 1°. In Fig. 9 we show the resulting surface pressure dis- 
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tribution where a comparison is also made with the results obtained by 
Jameson (Ref. 24) who used a steady relaxation procedure similar to that 
of Korn and Garabedian. The major discrepancy here in the two results is 
in the fulfillment of the shock jump condition and in the pressures just 
downstream of the shock. The desired condition downstream of the shock is 
indicated here by the dashed horizontal line, and it is seen that despite 
the fine mesh used, the result of Jameson does not capture the shock correct- 
ly. We shall discuss the treatment of shocks by the various steady relaxation 
procedures further' in the next section. 

The fourth example is an aft cambered profile of 3*9$ thickness ratio, at 
IS^ = O. 9 O and a = 0. The resulting surface pressure distribution is shown 
in Fig. 10. Here only the first order planar slab condition was used, so 
that, as in the first example, when a second order planar condition is em- 
ployed, one can expect a somewhat further upstream displacement of the shock. 

The shock capturing ability of the present procedure is further examined for 
an aft-cambered airfoil (Airfoil A) previously calculated and tested at high 
Reynolds number (30 x 10 /15 inch chord) at NAE (Ottawa) (See Ref. 15)* 

The earlier comparison of the experimental and calculated result at = O .85 
showed a remarkable agreement in the pressure distributions (see Fig. 11 ) 
everywhere except for the pressure along the aft portion of the lower surface, 
the location of the shock wave, and the pressures downstream of the shock, 
which understandably are the expected departures due to viscous effects. In 
the above comparison the angle of attack in the calculations was simply 
adjusted until the plateau pressures on both the upper and lower surface 
matched. Such a procedure could then be used to obtain an indication of 
the incidence correction due to wall effects for a highly porous perforated 
tunnel, if the viscous effects on the airfoil were negligible. Viscous 
effects, however, are not negligible for the aft-cambered peaky Airfoil A 
even at the high Reynolds numbers of 30 x 10^/chord. Leading edge separation 
and the change of the effective aft camber due to the not insignificant 
boundary layer displacement effects on the aft portions of both the upper 
and lower surfaces will significantly alter the plateau pressure levels. It 
is somewhat fortuitous that a simple angle of attack correction under such 
a circumstance enabled one to obtain the agreement shown in Fig. 10. With 
the present updated unsteady procedure a further calculation was carried out 
as the fifth example in which a preliminary attempt was made to determine 
how the upper aft portion of the airfoil (the viscous ramp of Ref. 20) must 
be modified to move the shock to the experimentally observed location and 
to obtain the empirically measured pressures aft of the shock. The procedure 
was modified such that the experimentally measured pressures aft of the shock 
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on the upper surface could be prescribed, replacing the prescription of the 
airfoil shape. The results of the calculations would then yield the surface 
pressures where the shape was prescribed and the effective shape of the air- 
foil where the pressures were prescribed. Such a result could then be used 
to obtain the upper surface boundary layer displacement. In the initial 
calculation the angle of attack was taken at 1.4° of the earlier calculations, 
but it was found that the alteration of the upper aft surface by the viscous 
ramp decreased the effective aft camber, lowering both the level of the 
under-pressures on the upper surface and the overpressures on the lower 
surface. The angle of attack was then simply increased to recover a match 
of the upper surface pressures. In Fig. 12 we show the resulting pressure 
distributions and the shape of the viscous ramp in terms of the slope dis- 
tribution. The plus symbols on the aft upper surface denote the prescribed 
pressures. For this case the resulting shock at the surface represents a 
strong oblique shock where the downstream Mach number is subsonic. If we 
take the initial increment 5-1° as the wedge angle of the viscous ramp, then 
this is to be compared to 6.2° which would be required by the jump conditions 
to produce the shock pressure rise of the calculations for the calculated 
upstream Mach number of 1.28. This fair agreement is indeed encouraging 
though of course the use of a still finer mesh might have produced a more 
convincing result. There still remains however a difference in the lower 
surface pressures of Fig. 12 and experiments since the lower surface boundary 
layer displacement effects and the upper surface leading edge separation 
effects have not been taken into account in the calculation. The effect 
of the lower surface displacement effects would tend to decrease further the 
effective aft camber, further lowering both the upper surface under pressures 
and the lower surface over pressures. On the other hand leading edge separa- 
tion would tend to decrease the suction levels over the upper surface as well 
as to produce the observed undulation in the supersonic pressures ahead of 
the shock. The net result must be that the angle of attack must be further 
changed to counterbalance these effects. If such effects could be calculated, 
then the difference of the resulting value of incidence and the geometric 
incidence in the experiment would then finally represent the incidence cor- 
rection due to wall effects if such effects have not been sufficiently 
aggravated to invalidate a simple incidence correction. Needless to say 
viscous effects can play a complicating role even at high Reynolds numbers 
(though admittedly Airfoil A is an extreme profile), making a direct comparison 
of experiments and calculations a major undertaking. 

In the above calculations the treatment of the "wake region" required an 
improvisation since the wake could not be simply treated as a rear stagnation 
streamline where a straightforward matching of the pressures and streamline 
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slopes was sufficient. Clearly across the wake a matching of the pressures 
is still a reasonably valid condition, but the streamline slopes clearly 
would not be equal. In the first attempt only a match of the pressures was 
used as a boundary condition on the sides of the slab representing the upper 
and lower surfaces of the wake. As probably should have- been expected the 
calculations resulted in an instability in the wake region, suggesting the 
need for a second crossing relation analogous to the matching of the streamline 
slopes in the case of the rear stagnation streamline. Such a relation no 
doubt could be obtained by simultaneously calculating the wake flow, but 
in lieu of this, a jump condition in the streamline slopes across the wake 
was postulated as follows. First the slope of the centerline of the wake was 
prescribed to vary linearly from the average of the streamline slopes at the 
trailing edge to a zero value at a prescribed distance downstream of the 
trailing edge and then taken to be zero thereon downstream. The slope of 
the upper surface of the wake (actually the displacement surface) relative 
to the above meanline slope was then postulated to be positive and just the 
negative of the lower surface slope, again taken relative to the meanline. 

For high subsonic flows it would not be unreasonable to assume that the flow 
upstream of the trailing edge was relatively insensitive to the preciseness 
of the above postulated wake, though in the present calculation this in- 
sensitiveness was not verified. 

Lastly as a further exercise in prescribing mixed conditions along the air- 
foil surface (either streamline slopes or pressures) we shall next investigate 
the possibility of evolving shockless flows or flows with weakened shocks by 
prescribing for example a continuous pressure distribution along the aft upper 
surface. Before we proceed to this task we note first some properties of 
shockless airfoil flows that severely limit what may be accomplished in this 
manner. First it is well known that along the sonic line (which is an isobar) 
the Mach waves "arriving" at the sonic line can only be expansion waves, while 
those "leaving" the sonic line must be compression waves. This then requires 
that only compression waves can impinge onto the airfoil surface, while only 
expansion waves can leave the surface when the flow is shockless. This 
requirement then mandates that the supersonic segment of the profile must 
have sufficient convexity to produce sufficient expansion effects to offset 
the effects of the incoming compression waves. Clearly, having any finite 
segment of zero curvature where the flow is supersonic would violate the 
above necessary condition for a shockless flow. Consider now the case of 
Airfoil A treated earlier. This airfoil has an upper surface having essen- 
tially zero curvature over much of the region where reduced convexity will be 
required, so it does nob represent a promising candidate airfoil to render 
shockless. On the other hand it does represent an effective airfoil, and it 



will be an interesting exercise to prescribe a continuous aft upper surface 
pressure distribution to see to what extent the shock wave can be weakened 
(or shortened) . An approximate guide as to how far upstream the upper 
surface modification should be incorporated can be obtained by noting the 
approximate domain of dependence of the portions of the shock that are to be 
weakened. Thus in the case of Airfoil A at M = O .85 and a = 1.4° we shall 

00 

prescribe an arbitrarily chosen continuous pressure distribution on the upper 
surface aft of the 30$ chord station. The resulting pressure distribution 
and the shape of the resulting humping is shown in Fig. 13 . Here the 
prescribed pressures are shown by the plus symbols. Shock weakening must be 
carried out only to an extent that the resultant drag in a real flow be- 
comes just tolerable, since lift will be lost in the process of adding com- 
pression waves to the upper surface. In Fig. l4 we next, show the resulting 
isobar pattern to illustrate the extent to which the shock has been weakened 
and shortened. Elimination of the stronger portions of the shock in the 
neighborhood of the surface not only reduces the entropy losses, but it also 
reduces significantly the adverse pressure gradients that would confront a 
boundary layer. In the above process a blunt base will invariably arise. 

Some bluntness of course is necessary to permit the removal of the boundary 
layer displacement effects, but in general further iterations in the shaping 
are required to evolve an aft upper surface configuration which does not pro- 
duce adverse viscous effects. 

Upper surface contouring to weaken the shock is not new, having been proposed 
earlier in Ref. 3 in the form of humped profiles. ' More recently Steger and 
Klineberg (Ref. 25 ) used the small disturbance steady relaxation procedure 
to calculate a number of interesting examples where also the surface pres- 
sures were prescribed either partially or entirely over the profile. Un- 
avoidably the use of the small disturbance procedures can lead to numerical 
difficulties at the leading edge, especially in a lifting case, which seldom 
can be resolved in a satisfactory way. Several of the examples calculated by 
Steger and Klineberg were perhaps unacceptably compromised in the process of 
avoiding the leading edge problem by their assuming that the flow over and 
upstream of the nose remained invariant even when aft profile modifications 
produced a significant change in the lift. 

In summary the results of the present section are gratifying in many respects 
and have yielded results that are reasonably viable. However it would be 
clearly premature to claim them to be "exact" inviscid results, since in no 
case was an assessment made of the effects of the remaining truncation errors. 
Particularly care must be exercised in the treatment of the shock where 
accuracy degrading diffusive damping must be employed to critically damp the 
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Figure 14. Isobars for Airfoil A With a Shock Attenuator Hump 
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over and under shoots about the shock that are characteristic with second 
order difference schemes. The Zierep post shock pressure cusp has features 
of an overshoot, so that it is essential that a sufficiently fine mesh be 
employed about the shock to capture the shock automatically and correctly 
when the oscillations are critically damped. Calculations on the whole are 
gradually becoming more routine, but despite the numerous examples computed 
to date, somehow annoying and time consuming surprises still manage to appear. 
Clearly we are still in the learning process. 

Finally with regards to the computing time required, typically 300 seconds 
were required for a given case on the CDC 7600 computer using 4107 mesh points 
and 1500 time planes starting with an incompressible flow. Such times are an 
improvement over the original unsteady procedure and are tolerable for use in 
planar calculations, but clearly they are still intolerable for practical 
three-dimensional examples. 


6. CONCLUDING REMARKS 

In the reported effort the flow was considered directly in the physical plane, 
since an important prerequisite for the finite difference analogue was its 
suitability for use in a procedure to incorporate viscous effects using a 
boundary layer concept, and for extension to steady three dimensional flows. 
Thus the use of the slab boundary and a quasi-planar condition along the 
parallel sides of the slab, introduced primarily to simplify the problem, are 
of particular convenience in applying profile modifications due to boundary 
layer displacement effects and for extension to three dimensional configurations. 

The problem here was formulated as an unsteady (hyperbolic) problem with the 
equations expressed in appropriate conservation form in the subdomains con- 
taining shocks, to conform to the guidelines suggested by Lax that insure 
the proper jump conditions across the shock discontinuities. In contrast to 
this, in the steady formulation where the equations are of mixed elliptic - 
hyperbolic character, and additionally in some cases even of second order, the 
concept of a weak solution has not been established, and there simply is no 
guide, comparable to Lax's contributions, to pose the problem properly to 
insure proper shock jump conditions. Thus, for example, it is not clear 
whether, in the steady case, different shock jump conditions would result 
depending upon the form of the equations used as in the unsteady hyperbolic 
case. Here in the latter case as noted by Lax (Ref. 4) and Courant and 
Hilbert (Ref. 10), it may be recalled that even if the continuity and momentum 
equations are used in the identical conservation form, different shock con- 
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ditions are obtained if the fourth equation is expressed, not as the con- 
servation of total enthalpy, but of entropy. 

Quite apart from being able to define a weak solution, there is still the 
task of obtaining the thus -defined weak solution by a numerical procedure. 

Here again we are indebted to Lax who has shown that, at least for an initial 
value problem, the proper weak solution can be obtained by solving a properly 
posed finite difference problem. The results of the present calculations 
indicate that Lax's conclusion can be extrapolated to the more complex airfoil 
problem, where the proper shock jump conditions and Zierep's post-shock ex- 
pansion (we shall discuss this later) can be obtained when a sufficiently 
fine mesh about the shock is used. 

To review the situation for the treatment of the shock, in Fig. 15 we have 
plotted the values of Mg* (based upon a*) at the airfoil surface downstream 
of the shock as a function of the values upstream for both the unsteady, 
calculations (Ref. 3 and the present results) as well as the steady relaxation 
calculations (Refs. 6, 7> 17 > 18, 24). Also shown are the curves reflecting 
the jump conditions Implied by .the unsteady weak solutions, as well as those 
from the Rankine-Hugoniot conditions consistently approximated. Consider 

first the unsteady results. The results of the present calculations as noted 
earlier do indeed fulfill the expected jump conditions, while the earlier 
results of Ref. 3 exhibit a scatter, and in most cases the values of M 2 are 
somewhat below the correct values. These latter discrepancies can be elimina- 
ted by simply decreasing the excessive diffusive damping which has caused a 
round off of the Zierep cusped pressure bucket as was demonstrated in an ex- 
ample of Section 5* (See also Fig. 15 ) . On the other hand the steady relaxa- 
tion results have a significantly greater scatter, and the implied mean curves 
fall considerably below the corresponding Rankine-Hugoniot values. Most 
probably the scatter is due in great part to an inadequate mesh spacing, since 
the use of the potential rather than the velocities as the dependent variable 
would make more stringent the requirements on the mesh size. However the 
cause of the large departures from the desired shock jump conditions remaining 
after the scatter is eliminated is uncertain at present. Though such a de- 
parture could be caused either by the excessive damping effects rounding off 
the Zierep cusped pressure bucket as in the unsteady case described earlier, 
or possibly to the form of the equations implying jump conditions other than 
the desired ones, it is also possible that the cause may be due to an inadequacy 
in the numerical procedure itself where only a dichotomy of difference schemes 
is employed to represent a continuous variation of the "domain of dependence. " 

In other words it may be conjectured that the flow is "shocked" not by the 
physics, but by the abrupt change in the difference scheme. Under these 
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conditions it would be unreasonable to expect the correct jump conditions. 
Clearly the picture would be clarified by carrying out additional calculations 
with refined mesh schemes. Thus at the present stage one must hold in abeyance 
the shock- capturing ability of the steady relaxation procedures. If a refined 
mesh does not remedy the situation* one should seriously consider a modification 
of the procedure to fit in a shock discontinuity directly, imposing the shock 
jump conditions, as was done earlier by Emmons (Refs. l6 and 17)? whose results 
(Fig. 7, Ref. 17) as shown in Pig .15 do indeed fulfill the required jump 
conditions. With such a modification the resulting steady relaxation proce- 
dures, especially those where the airfoil is mapped into a circle, would be 
difficult to improve upon from the point of view of both accuracy and computor 
costs. 

In summary with regard to the treatment of the shocks, the unsteady finite 
difference analogue . as used is based upon a reasonably sound mathematical 
background, and the results show that the proper shock jump conditions are 
fulfilled, and hence the correct shock location should then follow. On the 
other hand the steady relaxation procedures must await further effort before 
the treatment of shocks can be credibly accepted. Here, in the absence of a 
reliable guide for the correct formulation of the numerical analogue to assure 
the proper capturing of the shock, one should seek instead means to fit 
explicitly a shock discontinuity into the flow using the proper jump con- 
ditions across the shock. 


The discussions above regarding the surface shock pressure rise pertain to 
inviscid flow. These inviscid results are seldom observed in technically 
important experiments due to the interaction of the shock wave with the boundary 
layer, which in cases of interest is fully turbulent upstream of the shock. 

As the consequence of the interaction either an abrupt thickening or a separa- 
tion of the boundary layer takes place at the foot of the shock, and the 
resulting displacement effect of the boundary layer changes the shock at the 
surface from a normal shock in the inviscid case, to an oblique shock leading 
to a decreased pressure rise across the shock. Since the location of the shock 
is in part dependent upon the shock pressure rise, as well as the nature of 
the subsonic pressure recovery downstream of the shock, it is of paramount 
importance to develop a procedure to treat the above viscous interactions. 

In Ref. 20 a procedure has been suggested not only to evolve a fully determined 
system of turbulent boundary layer equations in the Von Karman integral format 
to treat a separated boundary layer, but also to solve the strong viscid-inviscid 
flow coupling. Here the artifice of a wedge-shaped viscous ramp has been 
introduced to simulate the boundary layer displacement effects. 
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In deter mining the surface shock pressure rise from experiments, in particular 
where there is negligible boundary layer thickening downstream of the shock, 
one must keep in mind the possible presence of an abrupt post shock expansion 
measured by Ackeret, Feldmann, and Rott (Ref. 21) and deduced by Emmons (Ref. 
16) . If the pressure taps on the model are insufficiently spaced, then such an 
expansion can be overlooked, and an erroneous pressure rise can be inferred 
by an improper interpolation of the measured pressures. 

The behavior of the flow at the foot of the shock has been analyzed by Zierep 
(Ref. 22) and by Ferrari (Ref. 23) by a local analysis about the foot of the 
shock. For convex surfaces they found that the shock curvature was logarithm- 
ically infinite at the surface, the shock changing from a normal shock at the 
surface to an oblique shock at an infinitesimal distance off of the surface, 
and that an abrupt expansion occurred oust downstream of the shock (in fact 
starting with an infinite gradient). The uniqueness of such local singular 
solutions, as pointed out by Ferrari, is difficult to establish, but the 
measurements of Ackeret as well as some existing calculations do indeed ex- 
hibit the principal features deduced by Zierep. On the other hand the local 
analysis does not of course establish the scale of the cusped pressure bucket, 
both in terms of its depth as well as its width, so that if these scales are 
sufficiently small, then it is conceivable that a finite difference result 
can satisfactorily capture the shock without exhibiting a pressure bucket. 

When a significant pressure bucket does arise, the extent that it can be 
captured by a finite difference calculation would depend upon the mesh size 
used, since the truncation errors lead to dissipative effects which tend to 
round off the cusp. In the case that the airfoil surface at the base of the 
shock has zero curvature or a concave curvature as in the viscous case with 
a viscous wedge such a post-shock expansion would not necessarily arise. 
(Clearly such a post-shock expansion does not arise in the case of a detached 
shock ahead of a symmetric wedge at zero incidence at slightly supersonic M^) . 

t 

In the inviscid steady relaxation calculations erroneous shock pressure rises 
(see Fig. 15 ) have been obtained that are reasonably close to those that would 
be expected in experiments with the viscous effects present. It is highly 
tempting to then suggest that viscous effects need not be incorporated, 
especially since a viable viscous procedure has not yet been demonstrated. 

There are . two fallacies inherent in such a point of view. First of all it is 
not assured by any means that the proper shock location would be obtained even 
with the use of the pseudo-shock, since the reduced shock pressure rise would 
be in error off the surface, and it is not accompanied by the proper shock 
obliqueness at the surface. Perhaps of more importance is the fact that the 
boundary layer displacement effects, essential in determining the effective 
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airfoil camber, would not be accounted for .in such an approach; and it is well 
known that a change of the effective camber, especially in an aft-cambered air- 
foil, can affect the entire pressure distribution. On the other hand the 
results using the steady relaxation procedure with its pseudo-shock undoubt- 
edly yield surface pressure distributions in most cases that are closer to the 
experimental results than the inviscid results with a true shock. The steady 
inviscid relaxation procedures therefore would serve as a valuable predesign 
tool even if viable boundary layer procedures were eventually developed. 

In closing the authors wish to express their appreciation to Mr. Ralph 
Carmichael, who served as the NASA contract monitor, for his enthusiastic 
support and assistance during the present study, particularly during the post 
midnight sessions associated with time-sharing computors. 
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